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Abstract. The LIBOR market model is very popular for pricing inter- 
est rate derivatives, but is known to have several pitfalls. In addition, if 
the model is driven by a jump process, then the complexity of the drift 
term is growing exponentially fast (as a function of the tenor length). In 
this work, we consider a Levy-driven LIBOR model and aim at develop- 
ing accurate and efficient log-Levy approximations for the dynamics of 
the rates. The approximations are based on truncation of the drift term 
and Picard approximation of suitable processes. Numerical experiments 
for FRAs, caps, swaptions and sticky ratchet caps show that the approx- 
imations perform very well. In addition, we also consider the log-Levy 
approximation of annuities, which offers good approximations for high 
volatility regimes. 



1. Introduction 

The LIBOR market model (LMM) has become a standard model for the 
pricing of interest rate derivatives in recent years, because the evolution of 
discretely compounded, market-observable forward rates is modeled directly 
and not deduced from the evolution of unobservable factors, as is the case 
in short rate and forward rate (HJM) models. See |MSS97j . |BGM97| and 
[Jam97j for the seminal papers in LIBOR modeling. In addition, the lognor- 
mal LIBOR model provides a theoretical justification to the market practice 
of pricing caps according to Black's formula (cf. |Bla76j ) . However, despite 
its apparent popularity, the LIBOR market model has certain well-known 
pitfalls. 

An interest rate model is typically calibrated to the implied volatility 
surface from the cap market and the correlation structure of at-the-money 
swaptions. The implied volatility from caplets has a "smile" shape as a func- 
tion of strike, while its term structure is typically decreasing. The standard 
lognormal LMM cannot be calibrated adequately to the observed market 
data. Therefore, several extensions of the LMM have been proposed in the 
literature using jump-diffusions. Levy processes or general semimartingales 



as the driving motion (cf. e.g. |GK03j . EO05 , |Jam99j ). or incorporating 



stochastic volatility effects (cf. e.g. |ABRn5j . [WZOBj . |BMSn9p . 
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The dynamics of LIBOR models are typically not tractable under different 
forward measures, due to the random terms that enter the dynamics of 
LIBOR rates. In particular, if the driving process is a diffusion process or a 
general semimartingale, then the dynamics of LIBOR rates are not tractable 
even under their own forward measures. Consequently, even caplets cannot 
be priced exactly in "closed form" (meaning, e.g. by Fourier methods), let 
alone swaptions and other multi-LIBOR products. In order to calibrate the 
model, closed form solutions are necessary, and these are typically involving 
approximations. 

The standard approximation is the so-called "frozen drift" approxima- 
tion; it was first proposed by |BGM97] for the pricing of swaptions and 
has been used by several authors ever since. The frozen drift approximation 
typically leads to closed-form solutions for caplet pricing in realistic LI- 
BOR models, see |E6n5| and (BMSMj. Although some authors ((BDBOT], 
[DBSOlj and |Sch02j ) argue that freezing the drift is justified in the lognor- 
mal LMM, it is shown that it does not yield acceptable results for exotic 
derivatives and longer time horizons, see e.g. |KSS02j . Therefore, several al- 
ternative approximations have been developed in the literature. In one line 
of research, |KSS02j and |DG05j have derived lognormal approximations 
to the forward LIBOR dynamics (for deterministic volatility structures). 
Other authors have been using linear interpolations and predictor-corrector 
Monte Carlo methods to get a more accurate discretization of the drift term 
(cf. e.g. [HJJOlj and IG/OOQ . We refer the reader to [J508] and [CBMOGl 
Ch. 10] for a detailed overview of that literature, some new approximation 
schemes and numerical experiments. Although most of this literature focuses 
on the lognormal LMM, |GM03bj and |GM03aj ) have developed approxima- 
tion schemes for the pricing of caps and swaptions in jump-diffusion LIBOR 
market models, based on freezing the drift. 

In this article, we consider a LIBOR market model driven by a Levy 
process and aim at deriving efficient and more accurate log-Levy approx- 
imations (compared to the "frozen drift" approximation, for instance). As 
a main result, we develop log-Levy LIBOR approximations which may be 
represented as a deterministic drift term plus a stochastic integral of a de- 
terministic function with respect to a Levy process. In particular, in the 
context of Monte Carlo simulation the drift term can be computed out- 
side the Monte Carlo loop, while the stochastic integrals can be computed 
efficiently for each trajectory. In contrast, standard Euler stepping of the 
original LIBOR SDE involves, for each LIBOR trajectory, an accurate com- 
putation of a complex-structured random drift term at each Euler step and 
is therefore significantly more time-consuming. Theoretical investigations 
as well as numerical experiments show that the log-Levy approximations 
are both fast and accurate when the LIBOR volatilities are not too high, 
and thus provide an effective alternative to simulation methods based on 
standard Euler discretizations. Finally, as a generalization of [GBM06j . we 



^In a previous unpublished manuscript by the first and third author |PS10| the efliciency 
of the standard Euler approach was improved to some extend also, but there was still a 
costly random drift involved. 
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derive log- Levy approximations for annuity terms, which allow for pricing 
options in high volatility regimes. 

The article is structured as follows: in section [2] we review the Levy- 
driven LIBOR model, in section [3] we construct the log-Levy approximations 
to the model and in section 2] we provide some error estimates. Section [5] 
demonstrates numerically the effect of the approximations, while section [6] 
deals with an approximation of annuities. The final section provides some 
recommendations on the construction of multi-dimensional Levy LIBOR 
models, while the appendices collect various calculations. 



2. Levy LIBOR framework 

Let = To < Ti < • • • < T/v < Tjv+i = T* denote a discrete tenor 
structure where 6i = Tj+i — T,, i = 0, 1, . . . , A'', are the so called day-count 
fractions. For this tenor structure we consider an arbitrage free system of 
zero coupon bond processes Bi, i = 1,...,A^ + 1, on a filtered probability 
space (fi, J^, (J"t)o<t<T* 5 IP*), where := Pat+i is a numeraire measure 
connected with the terminal bond i?7v+i- From this bond system we may 
deduce a forward rate system, also called LIBOR rate system, defined by 

Ht) := y (4^7^ - ' < t < T,, 1 < i < TV. (2.1) 

Li is the annualized effective forward rate contracted at date t <Ti for the 
period [Tj,Tj_|_i]. |Jam99j derived a general representation for the LIBOR 
dynamics in a semimartingale frame work. In this article we consider a Levy 



LIBOR framework as constructed by EO05 ; see also |GK03j and [BSllj for 
jump-diffusion settings. 

Consider a standard Brownian motion W in W^, m < N, a bounded 
deterministic nonnegative scalar function a{s), s € [0,T*], and a random 
measure fj, on [0, T*] x with P^j-compensator F{s, dx)ds, where jj, and W 
are mutually independent. Let H = (-ff (i))o<t<T* be a time-inhomogeneous 
Levy process with canonical decomposition 



H{t) = J y^dVF(s) + J J x{fi{ds,dx)- F{s,dx)ds). (2.2) 

OR™ 

We denote by Jl the compensated random measure of the jumps of H, that is 
Jl{ds, dx) := fi{ds, dx) — F{s, dx)ds. In order to avoid truncation conventions 
we assume that F satisfies the (stronger than usual) integrability condition 



//(|N|A|rf)F(,,d.)d,<co. 



M'" 

We further assume that 



J J exp (n"'"j;)F(s, dx)ds < CO, (2-3) 



||x|i>l 



4 



A. PAPAPANTOLEON, J. SCHOENMAKERS, AND D. SKOVMAND 



for all < (1 + e)M, with M,e > constants. Thus, by construction, the 
process {H{t))Q<t<T, is a P^,-martingale. The cumulant generating function 
of H(t), t G [0, T^], is provided by 

t 

ln]E[e"^^W] = J K,{u)ds, (2.4) 



where 

= ^||n||2 + J (e"^^-l-'u"rx)F(s,dx). (2.5) 

Along with the Levy martingale (j2.2|) we introduce a set of bounded de- 
terministic vector-valued functions Aj(s) G M™, i = 1, . . . , N, usually called 
loading factors. In order to avoid local redundances we assume that the ma- 
trix [Ai, . . . , A7v](s) has full rank m for all s E [0, T*]. Moreover, we assume 
that ||Ai(s)|| < M, for ah i, and || Ei^iC*)!! < ^> for all s e [0,r*]. 

The Levy martingale and the set of loading factors then constitute an 
arbitrage free LIBOR system consistent with ()2.ip . whose dynamics under 
the terminal measure P^, are given by 



t t 



Li{t) = Li{0) exp J b,{s)ds + J Xj{s)dH{s) , (2.6) 

\0 / 

i = 1, . . . ,N, where the drift terms in the exponent are given by 
..^-^alA^r-Elfel-A^A, (2.T) 




1 + 6jL,. 



Xjx\ Fi;dx); 



for details see [EO05]. For notational convenience, we set Lj_(s) := Lj(s—) 
in (|2.7p . while the time variable is suppressed. 

Due to the drift term (j2.7p . a straightforward Monte Carlo simulation of 
(|2.6p would involve a numerical integration at each time step, since the ran- 
dom terms ijg'.^i. appear under the integral sign. In order to overcome this 
problem, we will re-express the drift in terms of random quotients multiplied 
with cumulants of the driving process. We have that 

j=i+l ■' ■' 



p=l i<ji<---<jp<N ■'^ -'1 ■'P ■'P 



P+l 

g=l 0<ri<—<rq<p 
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the derivation is deferred to Appendix [XJ for brevity. Here k denotes the 
part of the cumulant k stemming from the jumps of L, that is 

K,(u) = j (e"^^' - 1 - u^x)F{s,dx). (2.9) 

Therefore, we can now avoid the numerical integration when simulating LI- 
BOR rates. However, another problem becomes apparent in this representa- 
tion: the number of terms to be computed in (|2.8p grows exponentially fast 
as a function of the number of LIBOR rates N , namely it has order 0{2^). 

Remark 2.1. In a practically applicable model, the loading factors Aj may 
be decomposed as follows: 

A,(t) = c,(7(Ti-t)ei„^(t) GR'", 

m(t) := inf{i : Tj > t}, ||ei|| = 1, ejej = pij, l<i,j<N, 

for constants q > 0, some (e.g. parametric) scalar function g > 0, and a 
correlation structure (pij) which resembles the correlations between forward 
LIBORs observed in the market. For instance, (pij) may be obtained as a 
rank-m approximation of a suitably parameterized full rank-A^ correlation 
structure; see |Sch05j for details. Further, the scalar function a may be taken 
as a constant that controls the influence of the Wiener noise with respect to 
the jump noise. 

Remark 2.2. Using semi- analytic pricing methods based on Fourier trans- 
forms, the Levy-driven LIBOR model may be calibrated to caplet volatili- 
ties for different strikes and maturities in the spirit of [BSllj . |EK07j and 
IBE.TPllj . 

Remark 2.3. The Levy-driven LIBOR model is constructed under the ter- 
minal measure Pat+i in this paper, for definiteness. As an alternative, for 
products with shorter maturity for instance, one may consider for some 
< T/v+i, a Levy-driven LIBOR model for t < under the measure 
P^, with respect to the numeraire bond B^. Another possibility is to con- 
sider as numeraire the spot LIBOR rolling over account 

^ , . m(t)-l 

BM := 1, Bo{t) := J] (1 + 5,L,(T,)), 

m{t) := min{m : Tm > t}, < t < T/v+i, 

and the numeraire measure Po associated with it. If one prefers to work in 
one of these other measures, the drift term (|2.7p has to be modified in the 
following way: for the Libor model in the measure P^, replace in (|2.7p . if 

i < N, the sum - ^^^+1 and the product Iljlj+i ~ Sj^i+i and Hjii+i 
respectively, and if i > A^ , by Y1]-n ^^^'^ ^/TVj-n I'espectively. Likewise, for 
a LIBOR model in the measure Po, replace in (|2.7p — X^jLj+i by '^'j=m{t) 
and the product Hj^j+i by l/n}=m(t)- refer to Jamshidian (1999) for 
more details. The proper choice of a numeraire measure under which the 
Levy-driven LIBOR model is constructed may depend on the set of LIBORs 
involved in a particular (structured) product which has to be evaluated by 
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simulation. In principle, one should choose the measure in such a way that 
the respective sum and product in the drift (j2.7p involve as few terms as 
possible. 



3. Efficient and accurate log-Levy approximations 

The aim of this section is to derive efficient and accurate log-Levy approx- 
imations for the dynamics of the LIBOR rates under the terminal measure. 
This is based on an appropriate approximation of the drift term, cf. (|2.7|) . 
which has two pillars: 

(1) expansion and truncation of the drift term, 

(2) Picard approximation of suitably defined processes. 

We will first provide an overview of the approximation argument, and then 
present the full details in some particular cases. 

3.1. Outline of the method. Let us denote the log-LIBOR rates by Gj. 
They are defined via 

Gi{t) := log Li{t), 
and satisfy the integrated linear SDE, see ()2.6p . 

t t 

Gi{t) = am + j h{s)ds + j Xj{s)dH{s), (3.1) 



< t < Ti, I < i < N . The semimartingale characteristics of Gi are 

B' = [ bi{s)ds 
Jo 

C'= [ |Ai|2(s)a(s)ds (3.2) 
Jo 

[ [ lA{x)F\s,dx)ds= [ [ lA{Xjis)x)Fis,dx)ds, 

Jo Jr Jo iM™ 

where ^ E 5(M\{0}). 

Inspired by the lognormal approximation developed by |KSS02j in the 
context of the lognormal LIBOR market model, we will derive log-Levy 
approximations for the dynamics of Lj , or equivalently Levy approximations 
for the dynamics of Gj. The standard remedy for the numerical problems 
arising in LMMs is to "freeze the drift" , that is to replace the random terms 
in ()2.7p - or ()2.8p - by their deterministic initial values. In the present 
model, this obviously leads to a log-Levy approximation, which however is 
not accurate enough. 

The method for deriving efficient and accurate log-Levy approximations 
we propose can be summarized in the following steps: 

• consider the different product terms ^■'^ . . . '^^f'^i'' =: Xo, ^ 

in ([23]), where i 1 < ji < • • • < ip < iV; 

• define functions h : — )• M such that 



h{Gj^,. ■■ ,Gj^) - X. 
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apply Ito's formula to Xj^ . j^, which leads to an SDE of the form 
dX,,...,-^(s) = A,,.„,^(s,L(.))d5 + S,„..,^(s,L(s))TdVF(5) 

+ / Cj^,„j^{s,x,L{s))Jl{ds,dx), (3.3) 



with L = [Li, . . . ,Ln]; 

• use the first step of a Picard iteration to approximate Xji...jp by the 
Levy process 

t 

= ^nM^) + J ^n-.A''^^ m)ds (3.4) 



t t 
+ / Bj,...j,{s,L{0)ydWs + J I Cj,...j^{s,x,L{0))Ji{ds,dx); 

M™ 

• plug the Levy processes xjj^ into bi, cf. ()2.8p . which leads to a 
Levy approximation for bf, 

• finally, integrate by parts to deduce a Levy approximation for Gi of 
the form 

t t t 

Gi{t) ^ Gi{^,t) + j H(t,s)ds + j Q^{t,s)dW{s) + j I(t, s, x)/i(ds, dx), 



where H, O and I are deterministic, time-dependent functions. 

The main advantage of the above approximations is that they can be sim- 
ulated efficiently, as explained in section Moreover, their characteristic 
functions can be given in closed form. 

Remark 3.1. Note that the "frozen drift" approximation can be easily 
embedded in this scheme. It corresponds to using just the initial values 
Xj^...jp(0) instead of the Levy process X^^^ in (|3.4|) . 

3.2. Log-Levy approximation schemes. In the sequel, we are going to 
follow this recipe for deriving efficient and accurate log-Levy approximations, 
and present the full details of the method. However, we will first truncate 
the drift terms at the second order, in order to reduce the number of terms 
that need to be calculated. 

1. The first step is to expand and truncate the drift term at the second 
order; these computations have been deferred to Appendix lAl for brevity, see 
(|A.5p . We will approximate bi by where 

i+l<k<l<N ' « I ( i 

where 

9i = K{Xi), r]ij = K{\i + Xj) - K{Xi) - K{Xj) (3.6) 
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and 

Cikl = + Afc + A;) - K{Xi + Afc) - K,{Xi + Xl) 

- K(Afc + Xl) + k(A,) + K(Afc) + k(AO. (3.7) 

The number of terms to be calculated is thus reduced from 0{2^) to 0{N'^), 
while the error induced is 

bi = b1 + 0{N^6^\\L\f). (3.8) 

Therefore, the gain in computational time is significant, while the loss in 
accuracy is usually relatively small. The numerical examples verify this, see 
section ISTT] for more details. 

2. The second step is to approximate the random terms 

in (|3.5|) by a time-inhomogeneous Levy process. Define the functions 

X (Jj-e^ , / X 4e'''= die""' 

jyx) = , ' ^ and g{xk,xi; 



l + (5je^ ^ ^ 1 + (Jfce^* 1 + (5ze^' ' 



where 



(l + 5.ex)2 / (i + 5.e-)3 ■ 

The partial derivatives of g can be computed equally easily, and are denoted 

d d 52 , , 

9k = T. — g, gi = -K—g, gu = ^ — -^g, (3.10) 

OXk 0X1 OXkOXi 

and so forth. We obviously have that 

Zj(.t) = f{Gj{t)) and Yki{t)=g{Gk{t),Gi{t)). (3.11) 

The functions / and g are C^-differentiable, hence we can apply Ito's 
formula for semimartingales (cf. e.g. |JS031 Theorem 1.4.57]) to Zj and Y^i. 
Using (j3.ip we may derive (with time variable s suppressed or denoted by • 
in the integrands) 

dZ, = (^J (/(G, + Ajx) - /(G,) - /' (G,) Ajx) F(.,dx) (3.12) 

+ /' (Gj) b'; + i/" (Gj) I A,f ds + /' (G,) V^XjdW 
+ I (^f{G,^ + X]x) - f{Gj^)) (^(ds, dx) - F{;dx)ds) . 

R™ 

The derivation is given in Appendix [Bj Hence, we have that 
dZj{s) = Aj{s,L{s))ds + Bj{s,Lj{s))dW{s) 

+ j Gj{s,Lj{s),x){iJL{ds,dx)- F{-,dx)ds), (3.13) 



APPROXIMATIONS TO LEVY LIBOR MODELS 9 

with obvious definitions of tlie deterministic functions Aj , Bj , and Cj . Due 
to the drift term b'j, the function Aj depends on the whole LIBOR vector L 
rather than Lj only. 

Similarly, we have for Y^i that 

dYkiis) = Aki{s,L{s))ds + Bl{s,Lki{s))dW{s) 

+ j Cki{s,Lki{s),x) {fiids,dx) - F{-,dx)ds) , (3.14) 

where Aki, B^i, and Cki are deterministic functions; see Appendix [C] for all 
the details. Analogously to ()3.13p . A^i depends on the whole LIBOR vector 
L, while B/ii and Cki depend on Lk and Li only; this is denoted by Lki- 

3. The next step is to approximate Zj and 1^/ by suitable Levy processes. 
This approximation is based on a Picard iteration for the SDEs in (j3.13p 
and (|3.14|) . Regarding Z, the initial value of the Picard iteration is 

while the first order Picard iteration is provided by 

t t 

Zf\t) = Zj(0) + j Aj{s,L{0))ds + j Bj{s,Lj{0))dW{s) 



Cj{s,Lj{0),x){fi{ds,dx) - F{-,dx)ds) . (3.16) 

We can easily deduce that Z^^^ is a time-inhomogeneous Levy process, since 
the coefficients ^j(-,L(0)), Bj{-,Lj{0)), and Cj{-, Lj{0), ■) in (IXT6|) are de- 
terministic. Indeed, we have that 

Aj{s, L(0)) = /' (G,(0)) bf\s) + i/" (G,(0)) |A,f {s)a{s) 

+ I [f{Gj{0) + X]{s)x)-f{G,m-f'{G,{0))X]{s)x)F{.,dx), (3.17) 



where 



i+l<j<N ^^'^^'^^-^^^ 
i+l<k<l<N ' K W > I I \ J 



and 



Bj{s, L,(0)) = /' (G,(0)) vM^A,(s), (3.18) 

Cj{s,L,{0),x) = f{G,{0) + Xj{s)x) - fiGjiO))- (3-19) 
Analogously, the initial value of the Picard iteration for (j3.14p is 
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and the first order iteration is 

t t 
Yll\t)=Ykm+ j Aki{s,L{Qi))ds + j Bli{s,Lkm)dW{s) 



+ j J Ckiis,Lki{0),x){fi{ds,dx)-F{;dx)ds), (3.21) 



and we can again deduce that Y^p is an additive Levy process. 

4. The fourth step is to apply the Levy approximations of the random 
terms to ()3.5p . Let us denote by bi the resulting approximate drift term; we 
have that 

6'/«6,:=-0, - Yl E C^fczn?^. (3.22) 

i+l<j<N i+l<k<l<N 

Keeping in mind that bi will be integrated over time, we define 
t t 

Vij{s,t) = J r]ij{r)dr, and Viki{s,t) = j Cikl{'r)dr, 
s s 

which are obviously deterministic processes of finite variation. Now, for fixed 
t > 0, we can apply integration by parts, which yields 

t t 

m,{s)zf\s)ds^ V,,{Q,t)Zj{{)) + j Vij{s,t)Aj{s,Lmds 



t 

+ 



t 

+ j Vij{s,t) j Cj{s,Lj{0),x)Jl{ds,dx). 



Similarly for the other term we get 

t t 
J Qki{s)Y^\s)ds^Viki{0,t)Yki{0) + JViki{s,t)Aki{s,L{0))ds 



t 

+ jv,ki{s,t)Bl{s,Lkm)dW{s) (3.24) 



t 

+ j Viki{s,i) j Ckiis,Lki{0),x)Jl{ds,dx). 



5. Finally, collecting all the pieces together we can derive a Levy approx- 
imation for the log-LIBOR rates. The approximate log-LIBOR is denoted 



z 

J V^jis,t)Bj{s,Lj{0))dWis) (3.23) 
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by Gi and has the foUowing dynamics 

t t 

G^{t) = Gi(0) + j hi{s)ds + j \J{s)dH{s), (3.25) 



which using ([321, and (fCTjl leads to 

t t t 

= G,(0,t) + j Hi(t,s)ds + j Qj{t,s)<lW{s) + j s, x)/I(ds, dx), 



(3.26) 

where 

G,(0,t) :=G,(0)- ^i^(0,i)^i(0) 

«+i<i<JV 

i+l<k<l<N 

R^it,s):=-9i{s)- Yl Vijis,t)A,{s,L{0)) 
- Y V,kiis,t)Aki{s,L{0)), 

i+l<k<l<N 

eJit,s) := y^Xjis) - Y Vrj{s,t)Bjis,Lj{0)) 

i+l<j<N 

- Y V,ki{s,t)Bj^{s,Lkim 

i+l<k<l<N 

and 

li{t,s,x) := Xj{s)x - Y yij{s,t)Cj{s,Lj{0),x) 

i+^<j<N' 

- Y Vikiis,t)Ckiis,Lki{0),x). 

i+l<k<l<N 

Let us introduce the process Xf^(r), 0<r<t, defined by 

r r r 

(r):=Gi(0,r) +J Ri{t,s)ds+ J QJ {t, s)dWis) + J Ii(t, s, x)/I(ds, dx). 



Obviously, xf^(r), < r < t is a time-inhomogeneous Levy process whose 
characteristic function may be expressed by the Levy-Khintchine formula 
in terms of Hj, 0j and Ij in a straightforward manner. 

Remark 3.2. We will call the approximation in (|3.26|) the second order 
log-Levy approximation of the LIBOR rate. If we ignore the second order 
terms (i.e. those depending on and L;), we immediately arrive at the 
first order approximation. The numerical results in section [S] document the 
improvement from the first to the second order approximation. 
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Remark 3.3. If we restrict our model to the Brownian motion case, the ap- 
proximation in (j3.26p coincides with the "fuhy lognormal model" of |DG05j ; 
see also |KSS02] . 

Remark 3.4. Note that the approximation methods developed in the pre- 
vious sections do not depend crucially on the choice of the measure. If we 
work under the spot measure, cf. Remark 12.31 then the Picard approxima- 
tions can be carried out similarly. However, an additional approximation is 
required to represent the drift in terms of cumulants as in eq. (|2.8p (because 
of the l/rij terms). 

3.3. Efficient simulation of the log-Levy approximation. In this sec- 
tion, we outline how simulation of the Levy approximation 

t t t 

Gi{t) = GiiO,t) + jRi{t,s)ds + J eJ{t,s)dW{s) + Jli{t,s,x)Jl{ds,dx) 



(3.27) 

can be carried out in an effective way due to the fact that Gj(0,t) and the 
integrands in ()3.27p are explicitly known deterministic functions. 

(I) The terms Gi{0,t) and jQHj(t,s)ds are deterministic integrals which 
may be computed outside any Monte Carlo loop using some quadrature for- 
mula. 



(II) The Gaussian part 



t 

Q(t) := J eJ{t,s)dW{s) (3.28) 



may be computed either by usual Euler stepping, or even directly at some 
fixed time t if only the distribution of G{t) matters. In this respect, the distri- 
bution of any vector (^^(t), ...,Q^{t)) — for simulating a set of log-LIBORs 
(Gi^ (t), Gjj. (t))) — is Gaussian with explicitly known covariance struc- 
ture, and thus can be simulated straightforwardly. 

(Ill) Finally, consider the practically important case where the Levy mea- 
sure itself is time homogeneous, i.e. F{dx) = F{-,dx). After truncating 
this measure with respect to jumps with size smaller than some e > (if 
needed), simulation of a realization of the jump term in ()3.27p may effec- 
tively be carried out as follows. First sample on the interval (0, t) the num- 
ber Nf (of jump times) according to a Poisson distribution with intensity 
tF({||x|| > e}). Next distribute Nt jump points {si, satj} uniformly over 
the interval (0, t), and sample independently for each jump point si a jump 
xi, 1 < ^ < -/Vt from the probability measure 

F(dxn {||x|| > e}) 
F{m\>e}) ■ 
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Then a realization of the (compensated) jump term is obtained as 

m * 

'iiit) ■■=Y,^i{t,si^,xi)- / li{t,s,x)F{dx)ds, (3.29) 

'=1 ||x||>. 

where the deterministic integral term can be computed outside any Monte 
Carlo loop by standard methods. Note that a realization of the whole log- 
LIBOR vector {q^ {t), . . . ,(^^{t)) will be computed using the same set of 
jumps isi,xi), I = l,...,Nt. 

The main benefit from the log-Levy approximation as outlined above, is 
the fact that for the simulation of a log-LIBOR vector {Gi{t), ...,GN{t)), 
the computation of the terms in (|2.6|) via (j2.8|) or (j3.5|) based on each re- 
alization of the Brownian motion and the jump process on a fine enough 
time grid is not required. This is in clear contrast to the Euler (or predictor- 
corrector) discretization of (|2.6p and (j2.8p . It is obvious that in view of the 
complex structure of (j3.5p only, such a simulation would require the (accu- 
rate enough) construction of a whole log-LIBOR system {Gi{tj), Giy{tj)) 
for < ti < ■ ■ ■ < tn '■= t involving the evaluation of the function b" at each 
grid point tj. In contrast, simulation of the log-Levy LIBOR approximation 
only involves the evaluation of (|3.29p at the jump times and the relatively 
efficient simulation of the Wiener integral ()3.28p inside a Monte Carlo loop. 



4. Error estimates 

In this section, we will provide some error estimates for the log-Levy ap- 
proximations in order to offer a theoretical justification for the proposed 
approximations. The error estimates are rather qualitative in nature, how- 
ever they allow for useful conclusions. 

In view of (|3.25p we have for the pathwise error of the (log-)LIBOR ap- 
proximation. 



Li{t) 



U{t) 



< exp 



G,{t)-Gj{t) 



< exp 



bi{s) - bi{s) 



ds 



thus we need to study the difference \bi — bi\. Since the main contribution 
of this error is due to the first and second order term in (2.7), we consider 
instead (see ()3.5p ) 



(1) 



i%i+ 

i+l<k<l<N 



\Cikl\ 



Let us assume for simplicity that a(s) = 1, and that and are (dimen- 
sionless) constants such that 



, max |7yi |<Jf max sup ||Ai(t)||2 =: i^r,A 
max \Ciki\ < Kt max sup ||Aj(t)||9 =: KtX 

l<i<k<l<N ' l<i<A^O<t<T 



2 

max J 



2 

max" 
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We then have 

L^{t) 



log 



Li{t) 



< K„X 



V max 



max 



max 



«+i<i<Af , 



L2(P*) 



i+l<k<l<N 

For the term (/) we get from (|3.13p and (|3.16p 





t 

l\\Yil\s)-YMl^^^^^ds=:iI)Hn). 



Zf\s)-Z,{s) 



i2(P.) 



< / |A,-(n,L(0))-^,(tx,L(n))|^^(p^) 



+ { j E\\Bj{u,Lj{Q)) - Bj{u,Lj{u))\\ldu 



1/2 



1/2 



E {Cj{u, Lj{0),x) — Cj{u, Lj{u), x))^ F{u, dx)du 



In view of ()3.17p . (|3.18p and ()3.19p . let Ka, Kb, Kq be dimensionless Lip- 
schitz constants such that for all 1 < j < and < n < T*, 

\Aj{u,y) - Aj{u,y') \ < /^aALx h - > 



\\Bj{u, yj) - Bj{u, yj)\\^< KBXmax \yj - y'j 
j {Cj{u, yj,x) - Cj{u, y'j,x)f F{u, dx) < i^c^max \yj - y'j 



Then, using 



Z)"{s) - Z,{s) ^^^^^ < KaX^,, I \\m - L(n)||2,^^(p^) du 



+ {Kb + Kc) A^ax I J E |L,(0)) - Lj{u)\^ du 



1/2 



we obtain the estimate 

t 



{I)<Xi,,K^KA I ij ||L(0)-L(7x)||2,i^(p^)dn ds 

\0 / 
t I s 

+ Xl^Kr,{KB + Kc) I ^^mjc^lj E\L,m-L,{u)\''du\ ds, 



and a similar expression may be obtained for the second term (//). 
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On an intuitive level we may interpret the estimates (/) and (//) in the 
following way: if we roughly consider that (the approximate squared vari- 
ance) E \ Lj{0)) — Lj{u)\ ^ A^ax^' then for (/) we obtain 

t s 





t I s 



1/2 



udu 



ds 



^/2, 



K,,KAt'/^ + ^A^^K, {Kb + Kc) 



and a similar result for (//). Hence, for some dimensionless constants Ki 
and K2, 



log 



Li{t) 



<K,{Xl,^tf^ + K,{Xl,J)\ 



Concluding, the log-Levy LIBOR approximations are extremely good as 
long as A^ax^ is small enough but, may become poor as soon as this product 
grows very large. This issue is confirmed in our numerical experiments. 



5. Numerical illustrations 

Throughout this section, we will consider a simple example with a flat and 
constant volatility structure. Similarly zero coupon rates are generated from 
a flat term structure of interest rates: -6(0, Tj) = exp(— 0.04-Tj). We consider 
a tenor structure with 6 month increments (i.e. 5i = ^). As stated in the 
introduction, the Brownian motion case is already well studied; therefore 
we set a = 0, thus limiting ourselves to the case where H \s a pure jump 
Levy process. We consider two univariate specifications, for simplicity. The 
first is a tempered stable or CGMY process (cf. |CGMY02j and |MY08p 
with parameters M = G = 1?>, Y = 0.25 and C = 48.4201, resulting in a 
process with mean zero and variance 1 (at t = 1), infinite activity and finite 
variation. The CGMY process has cumulant generating function defined for 
all u G C with \Ru\ < min(G,M), 

«cGMY(n) = r(-y)G^ {{^-■^Y-^ + ^} 

.r(-.)M^{(i.^)^-i-f}. (5.1) 

The necessary conditions are then satisfied for term structures up to at least 
10 years of length because M = min(G, M), hence YA^l | Ai| < 12 < M. Ex- 
act simulation of the increments can be performed without approximation 
using the approach in |PT06j . This approach can be used when simulating 
from (j3.ip with or without drift expansions, but cannot be employed in the 
case of the log-Levy approximation in (j3.26p where jump sizes are trans- 
formed in a non-linear fashion. Instead we employ an approximation where 
we replace jumps smaller than e with their expectation which is zero since 
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the jumps are compensated. This means that jumps bigger than e follow 
a compound Poisson process which can be easily simulated using the so- 
called Rosinski rejection method (see |Ros01] and |AG07l p. 338]). We set 
the truncation point sufficiently low, at e = 10~^, thus making the variance 
of the truncated term J^^x^v{dx) = 3.11 x 10~^, which can be considered 
small enough to safely disregard. To be consistent, we employ this procedure 
everywhere we simulate from the CGMY process. 

The second specification is a compound Poisson process with normally 
distributed jump sizes — often referred to as the Merton model. The cumu- 
lant generating function for u G C is 

A (exp(/iu + a'^u^ 



K-Mertoniu) 



l-fLu). (5.2) 
1/A yielding a process with mean zero and 



We set A = 5, /2 = and a = y 
variance 1 (at t = 1), as before. 

In order to verify the validity of our approximations we consider linear, 
nonlinear and path-dependent payoffs; in particular, forward rate agree- 
ments (FRAs), caplets, swaptions and so-called sticky ratchet caplets. To 
price FRAs and caplets with strike K maturing at time Tj, we compute the 
following expectations: 

N 

FRAo = 6iBN+i{0)lEp,[ Y[ {l + 5iLi{Ti+i)){Li{Ti) - K)\, (5.3) 



l=i+l 
N 



5iSAr+i(0)lEp,[ II {l + 6iLi{Ti+i)){Li{Ti)-Ky 



l=i+l 



(5.4) 



Following |Klu051 pp. 78], we have that the price of a payer swaption with 
strike rate K, where the underlying swap starts at time Tj and matures at 
Tm {i < m < N) is given by 



k=i 



N 



where 



-1, 

1 + hK, 



l=k 



k = i, 

i + 1 < k < m 
k = m. 



(5.5) 



1, 



(5.6) 



Similarly, a sticky ratchet caplet, which is a path-dependent derivative, can 
priced by computing the following expectation: 

N 

Mo = 5i5^+i(0)]Ep.[ J] (l + 5iLz(r,+i))(i?,(r,))+], (5.7) 



1=1+1 



where 



Ri{t) = Li{t) - min{Li(ri), . . . , Li_i(T,_i)}, Vt G [Ti, T^] 



Note that sticky ratchet caplets are often embedded in mortgages as a pro- 
tection against interest rates moving above a historical minimum value. 



APPROXIMATIONS TO LEVY LIBOR MODELS 



17 



5.1. Performance of the drift expansion. As we have argued in section 
13.21 the truncation of the drift term in equation (j2.7p is necessary in order 
to build a model that is computationally tractable. This section illustrates 
the effect of this truncation using the standard Euler discretization of the 
actual dynamics, i.e. equations (|2.6p and (|2.8p . 

Due to the complexity of calculating the true drift we limit ourselves to 
setting = 10, corresponding to a 5 year term structure. Furthermore we 
consider volatility structures constant and flat at Aj = 0.2 and Aj = 0.6 re- 
spectively. We simulate 10000 paths and plot the absolute difference between 
the prices from the drift expansions and the price without expansion (i.e. 
the full drift in ()2.7p ) in Figures l5.ll and 15.21 Each Monte Carlo simulation 
is done using the same random shocks for each method, thus eliminating 
the Monte Carlo noise as an error source. The figures demonstrate that the 
effect of the truncation depends mostly on the level of volatility Aj and less 
in the choice of product to price or the driving process. Furthermore, we 

ERRORS ATM FRAs CGMY h.0.2 ERRORS ATM FRAs Merton i .0.2 
0.14| ^ ^ ^ ^ ^ 1 1 I ^ ^ ^ ^ ^ 1 1 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

Maturity Maturity 



Figure 5.1. Drift expansion: low volatility scenario. 

notice that for low volatility even the first order expansion can be considered 
adequate, since the maximum of the absolute error is smaller than 0.2 bp. 
Conversely, for the high volatility case, the second order expansion is nec- 
essary to get proper accuracy. However, going to the third order expansion 
or beyond appears to be unnecessary as there is no visible gain in accuracy 
(< 10~^ bp). Hence, in the next sections we will use the second order drift 
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ERRORS ATM FRAs CGMY 1.0.6 ERRORS ATM FRAs Merton 1.0.6 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 



Maturity Maturity 
ERRORS ATM CPLs CGMY 1..0.6 ERRORS ATM GPLs Merton 1.0.6 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

Maturity Maturity 



Figure 5.2. Drift expansion: high volatihty scenario. 

expansion as our benchmark case since any resulting error is small enough 
to be disregarded. 

In Table 15. H CPU times are shown when simulating 10000 paths on an 
Intel 17 PC running Matlab. Here we can see that highly significant speed- 
up is achieved when truncating the higher order drift terms, whereas the 
decrease in speed when taking higher order approximations into account is 
relatively negligible. The CGMY is slower than the Merton model due to 
the much higher jump intensity needed in its approximation. We conjecture 
that the efficiency can be improved using the methods of |KHT10j , but this 
lies outside the focus of this article. 

Full Drift 1st order 2nd order 3rd order 
Merton 358.5 3^95 448 479 

CGMY 471.9 16.29 16.59 16.74 

Table 5.1. CPU Times (sees) for 10000 paths 

Finally, to conclude the subsection we should also mention that pricing 
errors for swaptions and ratchet caplets(not shown here) are of similar order 
of magnitude as in case of caplets. 

5.2. Performance of the log-Levy approximations. Next we study the 
performance of the log-Levy approximations. We increase the number of 
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rates to the more realistic setting of = 20 and consider the pricing of 
FRAs, caplets, sticky ratchet caplets and swaptions. We consider swaptions 
on swap rates over the periods {Ti,Ti + 3) years. Since we have estabhshed 
that errors from the drift expansion can be disregarded, we consider as the 
benchmark case the second order drift expansion studied in the previous 
section. In Figures 15.31 and 15.41 we plot prices from the frozen drift, the 
first and second order log-Levy approximations of section [3l and include 
the annuity approximation of the following section for completeness (for 
the path- independent derivatives). We use both the Merton and the CGMY 
model. We can observe that the frozen drift is consistently beaten by both 
the 1st and 2nd order approximation in both models and for all four prod- 
ucts. The 1st and 2nd order log-Levy approximations have a quite similar 
performance suggesting that second order approximation may not be neces- 
sary. Note that other parameter values (higher /lower intensity for Merton 
and fatter tails/slower tail decay for CGMY) have also been studied and 
again the results are qualitatively the same. 



ATM FRAs Merton 



ATM FRAs CGMY 1=0.2 




Maturity 
ATM CPLs Merton X =0.2 



-benchmark 
■ log-levy 1 st order 
-levy 2ncJ order 
annuity approx. 
- frozen 



Maturity 
ATM CPLs CGMY X =0.2 




-benchmark 
iog-levy 1 st order 
-levy 2nd order 
annuity approx. 
-frozen 




Maturity 



Maturity 



Figure 5.3. Prices for the Merton and CGMY models. 



Concluding, the log-Levy approximations offer an alternative to the Euler 
(or predictor-corrector) discretization of the actual dynamics which can be 
simulated faster and yields almost as accurate options prices. 
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ATM SWPTs Merlon X=0.2 ATM SWPTs CGMY i=0.2 




01234567 01234567 



Maturity Maturity 



STICKY RATCHET CPLs Merlon 1=0.2 STICKY RATCHET CPLs CGMY 1=0.2 




123456789 10 123456789 10 



Maturity Maturity 



Figure 5.4. Prices for the Merton and CGMY models. 

6. Approximation of annuities 

In the lognormal LIBOR market model, it is well documented that prob- 
lems may occur for high volatilities due to a proportionally large Monte 
Carlo variance in the annuity term used for discounting under the terminal 
measure, see [BevlOj and [GBM06] . Motivated by this numerical problem, 
we will derive an approximation of the annuity term in the spirit of |GBM06| 
§10.13]. 

Let us define the annuity term 

N 

Mt)= n i^+^jHt)), (6.1) 
j=i+i 

and consider the vector of log-LIBOR rates G = [Gj+i, . . . ,Gn]. We define 
a function / : M^"* R such that 

N 

(xi+i, ...,xn)=x^ n + 

j=i+l 

The partial derivatives of / are provided by 

fk(.x) = ^f{x)= J] (l + 5,e"0<5fce^^ = /(x^ 
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for alH + 1 < < A^, while we obviously have that 

f{G{t))=Mt). (6.2) 
Applying Ito's formula to f{G), we have that 

N * 

= Mt) = MO) j MGis-))dG,is) 

N t 
j,k=i+l Q 

+ Y,{ A/(G(s)) - fjiG{s-))AG,{s) i . (6.3) 

s<t [ j=i+l J 

Noting that the annuity is a P^,-martingale, we will focus on the martingale 
parts of ()6.3p in the sequel. Using (jS.ip and the fact that H is also a P^,- 
martingale, we get that the martingale part of the first summand is 

N * 



^ TV 

The second summand is omitted, while the final summands yields that 

E|a/(g(.))- E /,(G(.-))AG,(. 



SMAW-A.(,-)X:^|MfzL^AG,(.) 

s<t I i=j+i 



AT 



s<i I i=i+l ^ ■' 



Ms) - Ms-) - A{s-) E il'^'rV , A,(g)4/^(d^.dx) 
M.^ \ J-«H-J- 



t 



Ms) - Ms-) - Ms-) E T^/rT^^^(«)^^^(«'d^)d«' 

(6.4) 

where the quantity Ai{s) in the last two integrals should be understood as 

N 

Ms)= n (l + '^iexp{G,(s-) + Aj(s)x}) . (6.5) 

j=i+i 
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Collecting all the pieces together, we have that the annuity Ai satisfies 
the following integrated SDE 

t 



A,{t) = AiiO) + J Ais-)K{s-)dH{s) 



+ j j {Ai{s) - Ai{s-) - Ai{s-)Ki{s-)} fJi{AsAx) 

R*" 

or, equivalently 

t 

Ai{t) = Ai{{)) + j Ai{s-)K{s-)dH{s) 



Ai{s) 



+ / / Ms-) 



Ai{s-) 



1 — Aj(s— )3; > /i(ds, dx), (6.6) 



where 



N 



Ms-) = E I 



. l + 5jLj{s-) 



Ms)- 



(6.7) 



The solution of the SDE ()6.6p is the stochastic exponential, thus we get that 

(t t 
j Ki[s-)m[s) -^J AjAi{s-)ds 


+ /■/■[ Ms 



R"^ 
t 



l}Jl{ds,dx) {6.i 
+ 1 I fi{ds,dx) I , 



Ai{s-) 

log ; Ms) \ Ms) 



R™ 



Ai{s-) j Ai{s-) 



where again Ai{s) should be understood as in (|6.5p . By freezing the ran- 
dom terms in the drifts and jump sizes in the above dynamics we get an 
alternative approximation for the annuity term. Note that the resulting ap- 
proximation is also a log-Levy approximation. 

We can now use this approximation to price caplets and swaptions, noting 
that their respective payoffs can be written in terms of annuities: 

"A,(T,+i: 



Co = 5Ar+i(0)]Ep, 

So = 5jv+i(0)]Ep, 



Am 

/ m 



A^i{Ti)-{l + 5iK)A{Ti^ 



k=i 



(6.9) 
(6.10) 



where the c^'s are defined in ()5.6p . A similar expression can be derived for 
the sticky ratchet caplet. 
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CPL Prices Merton T.-1.5 CPL Prices Merton T.-2.5 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

Level of Constant Volatility Level of Constant Volatility 



Figure 6.1. Caplet prices as a function of volatility {N = 20). 



6.1. Performance of the annuity approximation. In Figure 16. H the 
quality of the various approximations is studied for a number of at-the- 
money caplets as a function of the volatility. As before we set the number of 
rates to = 20, and simulate 50000 paths for each volatility level. The plot 
is for the Merton model while the results are similar for CGMY. Using that 
at-the-money call option prices are increasing and roughly linear functions of 
volatility (see for example |Wil98j , |BS94j and [BFW04j for the case of non- 
Gaussian distributions), we can observe that only the annuity approximation 
produces sensible option prices at all levels of volatility. Moreover, even the 
benchmark case fails when volatility grows beyond 30%, meaning that the 
Monte Carlo simulation has failed to converge. The frozen drift fails at even 
lower levels of volatility, while the log-Levy approximations fail at a higher 
level, similar to the benchmark case. The annuity approximation works for 
all (higher) levels and also, as we have seen in Figures [5 . 3 1 and [H7^ for the low 
levels. One should therefore be careful when the average (across maturity) 
at-the-money implied volatilities are above 30% which is indeed the case in 
the current market for USD denominated LIBOR caplets where volatilities 
range from roughly 80% in the short end to 25% in the long end (source: 
Bloomberg). 

Moreover, in Figure [6T2] we observe that this problem becomes significantly 
less severe when limiting the number of rates to 10 with 6i = 1 instead of 
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CPL Prices Merton T.=5 



benchmark 




log-levy 1 st order 
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Figure 6.2. Caplet prices as a function of volatility {N = 10). 



20 with 6i = 0.5. Needless to say, limiting the number of rates is rarely a 
possibility in practice. 

In order to intuitively understand why this approximation performs bet- 
ter in the high volatility case than the other methods (e.g. the standard 
Euler scheme or the log-Levy approximations), let us just concentrate on 
the lognormal case. We have from (|6.8p that 



iog Mt) 



N 

E 



Xj ■ VtM + deterministic terms, (6-11) 



where N denotes a standard normal random variate. On the other hand, 
from (|6.ip . we get that 

N 

log Ai{t) = ^ 5j exp [\j ■ y/ij\f + random terms^ , (6.12) 

j=i+i 

where actually the method of approximation will only affect the random 
terms. We can easily conclude from (|6.11|) and (|6.12|) that the variance of the 
annuity approximation is significantly lower that the variance of the stan- 
dard representation, which results in the faster convergence of the Monte 
Carlo method. Thus, the annuity log-Levy approximation should be inter- 
preted as a variance reduction technique for the LIBOR market model. 

7. Economically meaningful multi-dimensional Levy measures 

VIA subordination 

Next, we reflect on the properties the driving process should have for 
practical applications and provide some recommendations. In an economi- 
cally realistic Levy LIBOR model the very structure of the Levy measure 
is important. Since, from an economic point of view, any jump in the daily 
rate typically affects all segments of the yield curve, we require in our mod- 
eling that, at a jump time, all the LIBORs jump, not only the first or second 
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half of the LIBOR curve for example. Moreover, this requirement should be 
fulfihed regardless of the structure of the loading factors Aj ; the latter may 
be inferred from some calibration procedure for instance. A natural way to 
meet this property is to take Levy measures which are absolutely continuous. 
In a jump-diffusion setting this can be easily established by taking as Levy 
measure the product of one dimensional absolutely continuous probability 
measures pi, i.e. 

iy{dx) = pi{dxi) ■ ■ ■ pm{dxm); (7.1) 



see [BSllj . In this paper we consider LIBOR models based on Levy processes 
with possibly infinite activity, thus having available flexible and realistic LI- 
BOR models possibly without Wiener part (i.e. a = 0). However, Levy mea- 
sures of infinite activity cannot be obtained by simply taking the product 
of a set of one-dimensional Levy measures of infinite activity. Nonetheless, 
we seek for absolutely continuous infinite activity Levy measures such that 
the entailed jump processes maintain certain (weak) independence proper- 
ties. Such measures may be constructed by Brownian subordination (see e.g. 
[UT04]) as outhned below. 

Let Vl^ be a Wiener process on R"*. The characteristic function of W{t) 
is given by 



E 



JzW{t) 



iW4 



z e 



We now consider a subordinator {St)t>o on M+, with Levy triplet (0,0, p), 
and with Laplace exponent H, i.e. 

E [e"^*] = e*^(") := exp J (e"" - 1) p(ds)) , n < 0. 

(0,oo) 

Then the m-dimensional process Y defined by 

Y{t) := W{St) 

has characteristic function 

E [e^^'^W] = E [e [e*^''^(^')|5Jl = E [e^'^^^)] = e*=(*(^)) 



exp 



exp 



^s'${z) 



(0,oo) 



t I 2 

(0,00) 



II^IF 



lU(ds) 



l)p{ds) 



--: exp [t^{z)] 



As a result, y is a pure jump martingale Levy process with Levy measure 
satisfying 



^"^11' - l) p{ds) = I (e'^'^ - 1 - iz'^x)i^^{dx). (7.2) 



$(z) = / I e 2 

(0,oo) 
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It is easily checked that 

oo 

v^{dx) = [ ^ e-^ll-llV(d^)dx, (7.3) 







{V2^ 

which is a measure with absolutely continuous support. 

Example 7.1. Let {St)t>o be the inverse Gaussian subordinator with 

p{ds) = ^l^^^,yds, and E [e"^'] = ,-2ctV^{Vx^-^) , 
Then, ()7.2p is known explicitly as 



Cj2 



= -2cV^ W A + — \\zf - VX 



e.g. see |CT04j . 



Example 7.2. Let (S't)f>o be a Levy subordinator with the following prop- 
erties: 



/9(dt) = Ce-4^Z)_y(G)l|i>o}dt, 



E{u) = 2C7r(-y) 



(G2 - 2n)^/2 cos (^y arctan (^^)) " 



where D is the parabolic cylinder function. Then, (j7.2p is known explicitly 
as the Levy exponent of the CGMY process, cf. (jS.ip . with G = M; see 
|MY08| . 

Remark 7.3. By taking in ([2^2]) F{s, dx) := z^^(dx) with given by ([7131) . 
the jump-part of ()2.2p is represented by the process y constructed above. 
It is easy to see that Y has uncorrelated components, although they are 
generally not independent. Indeed, Y{t) has mean zero and we have that 



E 



y(fc)(t)y(0(i)l = E [e \w'-''\St)W'-^\S, 



0, 1< k < l< m. 



Thus in contrast to the jump-diffusion situation in [BSllj where all compo- 
nents jump at the same time independently, here the components of Y still 
jump at the same time but in an uncorrelated rather than in an independent 
way. 

8. Concluding summary 

We have presented a tractable numerical approach to simulate trajecto- 
ries of a general Levy LIBOR model in an efficient way. By this method we 
construct efficient approximations to the computationally demanding drift 
term in the Levy LIBOR dynamics. We have shown that, due to these these 
approximations, we arrive at a significantly more accurate log-Levy approx- 
imation than the one obtained by the usual "frozen drift" approximation. 
The performance of the method is illustrated by several examples. The pre- 
sentation is embedded in a flexibly structured multi-factor Levy LIBOR 
model which allows for natural modeling of mutual LIBOR dependences 
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(via incorporating suitable correlation structures). As such the paper sup- 
ports practical implementations of Levy interest rate models that, until now, 
played mostly an academic role. 



Appendix A. Computation of the drift 

A.l. Full expansion in terms of cumulants. We will derive a repre- 
sentation for the integral term of the drift ()2.7|) which does not involve an 
integration over random terms. Let us denote the integral term by 



e^' 



^ I SjLj^ { e^^^ - 1 



j=i+i 



Observe that 

I 

3 = 1 l<j<l l<jl<j2<l 

+ X/ '^h'^h^h + ••• + wi-'-Wi 

^<jl<j2<j3<l 

I 

= 1 + Ysl{wi, ...,wi), 
p=i 

where Sp denotes the elementary symmetric polynomial of degree p in / 
variables, i.e. 

Sl{wi,...,wi) := J2 ^jV'^ip' 1<P<^- 
i<ii<-<ip<« 

Thus Bj may be rearranged as follows: 

N-i 

V^^-l-A7x)F(-,dx) + ^ / (e^>-l) X 

P=i 



^^"^ y 1 + Ci.+i- 1 + 5nLn- ' ' "^"^^ 

:= (/) + (//). 
Let us consider in (//) for p>l the term 
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With jo := I, we may write 



^J--l) (e^Ji-_i)...(e^l-_i) (A.l) 



(^1 _ e^'o^) (l - e^ii") •••(!- e^l") 



9=1 



(-If+l [1 + 



where 



p+i p+i 

•e ^""9 



g=l g=l 0<ri<---<rq<p 

/ T T \ 

g=l 0<ri<---<rq<p ^ ^- 



0{\\xf) 

p+1 

5=1 0<ri<---<r5<p 

Obviously, expression (jA.ip is of order 0(||x||^) for any p > 1, hence (!) it 
must hold 

P+i 

(?=1 0<ri<---<rg<p 

Therefore, we can deduce the following representation for the integral term 
e^7^-l-A7x)F(-,dx) 

N-i r- J ^ J P+1 

^ 1 + 5,-, L,i- ' " 1 + Li _ ^ 

p=l i<jl<-<jp<N ^ ^ > q=l 

A,,,, +-+A, 



X E / ( e^-^-^^-^^-^^ ^ - 1 - (A,.^ + . . . + A,g X ) F(,dx) 

0<ri<---<rq<p 



^(A.) + E E T 



^jiLji- ^jpLjp- 



p=i i<ii<...<i^<iv - + ^n^n- 1 + -^ip^ip- 
P+1 

J^(_l)P+.+i ^ ^(a,.^+ + (A.2) 

g=l 0<ri<---<rq<p 



A.2. First order expansion of ()A.2|) . Let us consider the first order ex- 
pansion of Mi ; we get 

2 

i<j<n ^ ^ q=l 0<ri<---<ro<l 



+ 0(l|if)- 
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Note that 

2 



= - E ^(^>i) + E ^(A>i + A>J 

0<ri<l 0<ri<r2<l 

= (Ajo ) - (Aji ) + K ( + ) . 

Thus we obtain the following expression for the first order expansion of the 
integral term Mi 

i<j<N ^ ^ 

which leads to the following approximation for the drift term bi in ()2.7p 

b[ = K{X,)+ -Mi_(^(A, + A,)-K(A.)-K(A,)), (A.4) 

i<j<N ^ ^ 

taking also the terms stemming from the diffusion into account. 

A.S. Second order expansion of ()A.2p . Analogously, we can also derive 
a second order expansion of B j ; we get 

M, = K{Xi)+ V / 7 {K{Xi + Xj)-K{X,)-K{Xj)) 

i<j<N •> ^ 



+ E /'i'; faA. + A, + A,)-K(A. + A,) 

- 2{Xi + A;) - K{Xk + Xi) + K{Xi) + K(Afc) + k(A/)) 



which leads to the following second order expansion of hi in ([27 
6'/ = k(A.)+ E i + f.L. HA. + A,)-k(A,)-k(A,)) 

i<j<N ^ ^ 

+ E faA. + A., + AO-K(A. + A.) 

i+l<k<l<N 

- K{Xi + Xl) - K{Xk + Xl) + K{Xi) + K{Xk) + k{Xi)) . (A.S) 



Appendix B. Derivation of (j3.12p 

Using the ltd formula for general semimartingales (cf. jJS031 Theorem 
1.4.57]) we have 

Zj = Zj{0) + J f\G,{s-))dG, + ^ I r(G,)d(G,^G,^) 



+ E ifiGM) - f (G,(s-)) - /' (G,(s-)) AG,{s)) , (B.l) 

0<s<- 
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where (G^, G^) denotes the quadratic variation of the continuous martingale 
part of Gj, that is 

d(Gj^, G'j){s) = |Ajf (s)a(s)ds. (B.2) 
The sum in (jB.ip . using ()3.2p . may be written as 

{f{Gj{s-) + AG,{s)) - f {G,{s-)) - r (G,(s-)) AG,(s)) (B.3) 



0<s<- 



R™ 



+ 


Moreover, 



(/(G,(s-) + \]x) - f (Gjis-)) - f {G,{s-)) Ajx) f,{ds, dx) 
f{Gj{s-) + Ajx) - / {Gj{s-)) - f {Gj{s-)) \]x) F{s, dx)ds 
fiGjis-) + Xjx) - f {G,is-)) - f (G,(s-)) Ajx) ^(ds, dx). 



J /'(G,(s-))dG, = J f\G,{s-))b,ds + J f'{Gjis-))V^XjdW 



+ / / f'{Gjis-))XjxJlids,dx). (B.4) 



R™- 



Finally, by plugging (IM . dES]) and (|R4l) into (ITO . (fXT2D fohows. 

Appendix C. Derivation of (j3.14p 

The computations are completely analogous to Appendix[Bl thus omitted 
for brevity. The coefficients of 1^/ in (|3.14p are 

Aki{s,L{s)) = 9^{Gk{s-),Gl{s-)Ms) 

i=k,l 

+ ^ j 9ijiGkis-),Gi{s-))Xiis)Xj{s)ais) 
+ j (^g{Gk{s-) + Xlx,Gi{s-) + X]x)-g{Gk{s),Gi{s)) 

R™ 

- E 9^{.Gk{s-),Gl{s-))XJx\F{sAx)^ (C.l) 

i=k,l ^ 

Bl{s,Lki{s)) = J2 9^{Gk{s-),Gl{s-)),Ms)XJ{s) (C.2) 



i=k,l 
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and 

Cki{s,Lki{s),x) = ^ gi{Gk{s-),Gi{s-))Xj{s)x 

i=k,l 

+ g{Gkis-) + Xjx, Glis-) + Xjx) - g{Gkis),Gi{s)) 
-^gi{Gkis-),Gi{s-))Xjx. (C.3) 

i=k,l 

Here Lki{s) := {Lk{s), Li{s)) and denotes that Bki and Gki depend on Lk 
and Li (via Gk and Gi). 
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